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We study the steady-state behavior of a driven non-equilibrium lattice gas of hard-core particles 
with next-nearest-neighbor interaction. We calculate the exact stationary distribution of the pe- 
riodic system and for a particular line in the phase diagram of the system with open boundaries 
where particles can enter and leave the system. For repulsive interactions the dynamics can be 
interpreted as a two-speed model for traffic flow. The exact stationary distribution of the peri- 
odic continuous-time system turns out to coincide with that of the asymmetric exclusion process 
(ASEP) with discrete-time parallel update. However, unlike in the (single-speed) ASEP, the exact 
flow diagram for the two-speed model resembles in some important features the flow diagram of real 
traffic. The stationary phase diagram of the open system obtained from Monte Carlo simulations 
can be understood in terms of a shock moving through the system and an overfeeding effect at 
the boundaries, thus confirming theoretical predictions of a recently developed general theory of 
boundary-induced phase transitions. In the case of attractive interaction we observe an unexpected 
reentrance transition due to boundary effects. 



PACS: 05.70.Ln, 64.60.Cn, 02.50.Ga 
I. INTRODUCTION 

One of the basic properties of many driven, interact- 
ing many-body systems is the occurrence of shocks. A 
shock in a system of classical flowing particles marks the 
sudden transition from a region of low density to a re- 
gion of high density. A well-known example for a shock 
is the beginning of a traffic jam on a motorway where 
incoming cars (almost freely flowing particles in the low 
density regime) have to slow down very quickly over a 
short distance and then form part of the (high-density) 
congested region. A remarkable feature of such shocks is 
their stability over long periods of time, i.e., they remain 
localized over distances comparable to the size of parti- 
cles. In some sense one may regard shocks as soliton-like 
collective excitations of the particle system. 

Lattice gas models have proven to be excellent sys- 
tems for the theoretical investigation of shocks and of 
the consequences of shocks for the collective behavior of 
the particle system. An interesting situation is the cou- 
pling of such a driven lattice gas system to external par- 
ticle reservoirs . It is intuitively clear that unlike in 
equilibrium systems, here the boundaries will play a deci- 
sive part in determining the bulk behavior of the system: 
Since the system is open at the boundaries, particles will 
flow in and out and the current will carry boundary ef- 
fects into the bulk. Among the fundamental questions 
to ask in such a set-up is the stationary (i.e. long-time) 
behavior of the system as a function of the boundary 
densities. Numerical observations and mean-field based 
arguments show that varying boundary densities leads to 
boundary- induced phase transitions To understand 
why this happens it is clearly necessary to get insight into 
the collective behavior of the lattice gas and to investi- 
gate the role of the shocks. 



The best-studied example is the one-dimensional 
asymmetric simple exclusion process (ASEP) [^-^, a 
hard-core lattice gas where each lattice site can be oc- 
cupied by at most one particle and particles hop stochas- 
tically with constant bias to vacant nearest-neighbor 
sites. For this model not only the structure and mo- 
tion of shocks is largely understood, but also the sta- 
tionary phase diagram of the open system 1^ and the 
exact particle-density profiles are explicitly known 
as function of the external reservoir densities. Based on 
the exact solution of the ASEP and of a related lattice 
gas model [|j the role of shocks for the stationary phase of 
more generic open driven diffusive systems was elucidated 
This has led to a theory of boundary-induced 
phase transitions for one-component driven particle sys- 
tems in one dimension |^ , reviewed below. Of particular 
interest are systems where the stationary current j (p) as 
a function of particle density p has a single maximum, 
an example being traffic flow on single-lane or multi-lane 
highways The theory predicts a phase diagram 

with a flrst-order (non-equilibrium) phase transition be- 
tween a low-density phase (LD) and a high density phase 
(HD) and a continuous (second-order) order transition 
from both phases to a maximal-current phase (MC), see 
Fig. 2 below. In the context of traffic flow the proper- 
ties of the second-order transition suggest more efficient 
control mechanisms for avoiding jams Jl^ . 

The strength of the theory lies in the small number 
and generality of its assumptions. Hence it is tempting 
to test its validity in real traffic flow. The openness of the 
system models in- and outflow of cars on a road between 
two junctions. Indeed, the main features of the predicted 
flrst-order transition have been observed using data col- 
lected on a German motorway ||lj,|l^ . The second-order 
transition has been confirmed by Monte-Carlo simula- 
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tions of a suitably modified Nagel-Schreckenberg model 
for traffic fiow, originally introduced only for periodic 
boundary conditions . Since to our knowledge the 

data necessary to test the existence of the second-order 
transition in real traffic are not available at present, inde- 
pendent investigations of other traffic flow models are re- 
quired to establish the presence of such a transition. We 
consider here a continuous-time exclusion model with ad- 
ditional short-range interaction on top of the pure hard- 
core repulsion of the usual ASEP. The model is designed 
to be exactly solvable like the ASEP (to some degree) on 
the one hand and to be somewhat more realistic in de- 
scribing real traffic than the ASEP on the other hand. 
However, given the wide range of experimental appli- 
cations of hard-core lattice gases which comprise phe- 
nomena as diverse as the kinetics of protein synthesis 
I [I, diffusion in thin channels |l9[| , or polymer reptation 
1 2G| , broader relevance of our model is to be expected. 
Indeed, in order to obtain a more complete picture we 
also analyse a similar model with attractive next-nearest- 
neighbour interaction. This model is not related to traffic 
flow, but describes driven hard-core particles with attrac- 
tive short-range interaction. To our knowledge, this is the 
first investigation of the steady-state selection of a driven 
diffusive system with this type of interaction. 

The paper is organized as follows. In Sec. ^ we com- 
ment on some of the requirements of traffic flow modeling 
and we describe our model. This model is a special case 
of the driven diffusive systems studied some while ago by 
Katz, Lebowitz and Spohn ||2l[]. We explain in which re- 
spect this model is more realistic for traffic flow than the 
standard ASEP which is a limiting case of our lattice gas. 
This phenomenological explanation is then conflrmed by 
the calculation of the exact flow diagram, i.e. the sta- 
tionary bulk current as a function of the density. In this 
section the emphasis is on bulk properties and hence we 
investigate the system with periodic boundary conditions 
in the thermodynamic limit. In Sec. [II we review some 
details of the theory of boundary-induced phase tran- 
sitions of Ref. ||l^ and we introduce the boundary dy- 
namics for modeling the coupling of a finite system to 
boundary reservoirs of constant density. We give an ex- 
act solution for the stationary distribution along the line 
in the phase diagram corresponding to equal boundary 
densities (with proof postponed to the appendices) and 
discuss a mean-field analysis of the full phase diagram. In 
Sec. I^we present a mean-field analysis of the full phase 
diagram and discuss Monte Carlo data for the phase tran- 
sition lines while in Sec. ^ we perform a similar analy- 
sis for the model with attractive interaction. Finally we 
summarize our flndings and present some concluding re- 
marks (Sec. VI). Technical details of the derivation of 
exact results are presented in the appendices. 



II. ASEP WITH NEXT-NEAREST-NEIGHBOR 
INTERACTION - SOME COMMENTS ON 
MODELLING TRAFFIC FLOW 

The ASEP is a stochastic lattice gas of hard-core par- 
ticles with biased particle hopping. Particles hop with 
exponential waiting-time distribution with parameter 1 
to their nearest-neighbor site to the right, provided this 
site is empty. If the site is occupied, the hopping at- 
tempt is rejected. Symbolically one may represent these 
stochastic dynamics as follows 



A$ 0yl with rate 1. 



(1) 



Here A represent a particle and represents the vacant 
neighboring site. Even though this stochastic process 
is too simple to be a realistic model for traffic flow, 
some qualitative features of real traffic |2^j2^ can al- 
ready be seen: Shocks exist and the stationary current 
j{p) = p(I — p) as a function of particle density p has a 
single maximum. An apparently unrealistic feature of the 
particle distribution is the absence of correlations in the 
steady state which are seen in more sophisticated traffic 
flow models |17|Tll. An unrealistic feature of the current- 



density relation (known in traffic engineering as flow di- 
agram) is the reflection symmetry w.r.t. the maximal- 
current density p* = 1/2 and its rounded shape close to 
the maximum. 

The basic mechanisms which determine traffic flow ap- 
pear to be flrstly the competition between the desire 
to reach an optimal speed while keeping a (velocity- 
dependent) safety distance. When the traffic density be- 
comes sufficiently large this distance cannot be kept any- 
more and the speed has to be reduced. As a consequence 
the current drops at some density p* . Secondly, there 
is a certain amount of randomness due to variations in 
individual drivers behavior. This "noise" necessitates a 
statistical description of traffic flow phenomena. 

Various one-dimensional lattice gas models which in- 
corporate these mechanisms in different manners have 
been proposed The following is part of the picture 
that emerges: 

(i) The existence of a shock is generic and appears to 
be the consequence of the non-linear current-density re- 
lation g. 

(ii) The symmetric shape of the current-density relation 
results from particle-hole symmetry and is a feature of 
models in which cars move with constant probability or 
rate, independently of the environment beyond the near- 
est neighbor site to which they move. This is an un- 
realistic assumption since clearly car drivers slow down 
when they see a slowly moving car already some distance 
ahead. They do not just perform an emergency stop 
when the car is immediately in front of them. Numerical 
and analytical results for models |l^jlj,|2^ which allow 
for a reduction of speed that depends on the occupa- 
tion of sites further ahead show an asymmetric current- 
density relation resembling the shape of the current- 
density relation of real traffic. In these models speed 
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is implemented by jumps over a variable number of lat- 
tice sites. 

(iii) The unrealistically round shape of the current- 
density relation at p* is specific for the ASEP. Deter- 
ministic exclusion processes with parallel update 25 - 3C|] 
also show a symmetric current-density relation with one 
maximum, but the derivative of the current is discontin- 
uous at the maximal-current density p*, in this respect 
resembling the shape of the current in real traffic |2^] 
and of more realistic traffic flow models Increas- 
ing the hopping probability in a discrete-time process 
towards deterministic hopping, leads to an increasingly 
sharp jump in the current derivative at p* p7) . Hence 
the rather broad exponential waiting-time distribution of 
the standard ASEP seems to be responsible for the round 
shape of the current-density relation at p* . In terms of 
the motion of a single particle these dynamics correspond 
to an overestimated single-particle diffusion coefficient. 

(iv) For parallel update, but not for sublattice parallel 
update, increasing the hopping probability strengthens 
antiferromagnetic particle correlations ||l8| ^ ^8|| , i.e., cars 
are less likely to be found close to each other than some 
distance apart. 

In order to further investigate this picture and to dis- 
entangle the various effects associated with a small dif- 
fusion coefficient and with speed-reduction respectively 
it would be interesting to study both ingredients sep- 
arately for models in which the stationary distribution 
can be calculated exactly. A small single-particle diffu- 
sion coefficient can be implemented in discrete-time mod- 
els choosing the hopping probability close to one (low 
noise) . In the single-speed ASEP this has been shown to 
lead to the (almost) discontinuous behavior of the cur- 
rent derivative discussed above. The numerical results on 
the effect of a small diffusion coefficient on correlations 
are inconclusive. In our toy model we keep the exponen- 
tial waiting-time distribution (large diffusion coefficient), 
but introduce a next-nearest-neighbor interaction which 
in the repulsive case models slowing down of a car if the 
next-nearest-neighbor site is occupied as well. A particle 
hops to the right with rate r if the next-nearest-neighbor 
site is empty and with rate q if it is occupied: 
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with rate r 
with rate q. 



(2) 
(3) 



The condition q < r models slowing-down, in the limit- 
ing case r = q one recovers the usual ASEP. For q > r 
this model has not an interpretation as traffic model, but 
may be regarded as describing hard-core particles with 
attractive short-range interaction which are driven by an 
external field. 

This model is a special case in the class of driven dif- 
fusive systems investigated in Ref. On a ring with 
N sites with periodic boundary conditions the stationary 
distribution turns out to be given by the equilibrium dis- 
tribution of the one-dimensional Ising model. Each state 
of the system is defined by the set of occupation num- 



bers n = {rii, . . . ,nj^} with rii = 0, 1. The stationary 
probability of finding a state n is given by 



(4) 



Here is the partition function and the "chemical po- 
tential" h parametrizes the conserved bulk density p. 
This grand-canonical distribution is a non-equilibrium 
stationary state, i.e., it is invariant under the stochastic 
time evolution, but it does not satisfy detailed balance 
with respect to the dynamics. It is interesting to notice 
that correlations are non-vanishing and, in the repulsive 
case q < r which corresponds to speed reduction, be- 
come antiferromagnetic. In fact, the stationary state is 
identical to that of the discrete-time ASEP with parallel 
update for hopping probability p ~ I ~ q/r. 

According to the dynamics described above the local 
density satisfies the continuity equation 



dt 

with the current 



{it ) ("j(l - ni+i)[qni+2 + r{l - ?^^+2)] ) 



(5) 



(6) 



between sites (i,i -I- 1). The stationary particle current 
j — (ji) is constant and is readily calculated using stan- 
dard transfer matrix techniques for the one-dimensional 
Ising model js^ (see Appendix A) . In the thermodynamic 
limit N oo one finds the exact current density relation 



J =rp 



v/l-4p(l-p)(l-g/r)-l 



2(1 - p){l - q/r) 
shown in (Fig. |l|) in the repulsive case. 

Repulsive case 



(7) 
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FIG. 1. Stationary current j as a function of the density p 
for r = 1, g = 0.1 (repulsive interaction). 
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III. OPEN BOUNDARIES 



A. Theory of boundary-induced phase transitions 



The analysis of the previous section yields the bulk 
properties of the lattice gas at a given density. Now 
we address the question of the steady-state selection of 
the system with open boundaries, i.e., we investigate the 
bulk density of an open system coupled at both ends to 
reservoirs of different fixed boundary densities. Trans- 
lated into traffic language, open boundaries correspond 
to traffic junctions at the ends of a road where particles 
(cars) enter or leave the road respectively with certain 
fixed attempt rates, a situation envisaged with different 
aims already in |Q for the ASEP with parallel update. 
As a result of the coupling to reservoirs a non-trivial sta- 
tionary density profile in the vicinity of the boundaries 
will emerge and the (spatially constant) bulk-density will 
be a function of the two reservoir densities. 

There are two distinct mechanisms at work. Firstly, 
because of the particle interaction, coupling of a semi- 
infinite system to a reservoir will generically lead to 
some discontinuous behavior of the stationary distribu- 
tion close to the boundary. The boundary represents an 
inhomogeneity of the system since the interaction of the 
particles with the fixed boundary leads to different dy- 
namics than that which results from the interaction of 
particles among themselves. This is a non-universal phe- 
nomenon which depends on the precise nature of the cou- 
pling mechanism and on the nature of the particle inter- 
action. For short-range hopping and systems with short- 
range stationary bulk correlations one expects the follow- 
ing picture: Coupling of a semi-infinite system at site 1 
to a reservoir of constant density pL will give rise to a 
non-universal boundary density profile starting at pL and 
approaching (on the scale of lattice units) some bulk den- 
sity p_ which is a non- universal function of (Fig. H). 
A similar picture holds for coupling of a semi-infinite sys- 
tem at the right boundary where particles flow out of the 
system into the reservoir. Here the bulk density p+ may 
change close to the boundary to the reservoir density pR. 
Superimposed on this non-universal boundary structure 
is a universal behavior which depends only the effective 
boundary densities p-,p+ close to, but not at the bound- 
aries. The theory of boundary induced phase transition 
describes the stationary phase diagram in terms of 
these effective boundary densities. 




I I I I I I I I I I I I I I I I I I I 
1 i 

FIG. 2. Coupling of a semi-infinite driven particle system 
to a reservoir of constant particle density pL- The (interpo- 
lated) density profile pi (i — 1,2, . . .) approaches some con- 
stant value p- after some finite distance from the boundary. 



We review only the principal ideas of this theory |^ 
which is based on an interplay of the collective velocity 



dj 
dp 

of the lattice gas and the shock velocity 

j+ - 3- 



P+ 



(8) 



(9) 



of a shock with limiting densities p+ and p_ and with 
limiting currents j-^- and j_ to the right and to the left 
respectively. 

The collective velocity is the velocity of the center of 
mass of a local perturbation in a homogeneous station- 
ary background (Fig. ^). It is positive for background 
density p < p*, but becomes negative for p > p*. In 
this case the perturbation creates a back-moving traffic 
jam, which leads to a negative center-of-mass velocity, 
even though all individual particles move with positive 
velocity. In terms of traffic flow one might think of such 
a perturbation as being a car which has just entered a 
major throughway from some side road. 

The shock velocity describes the motion of the shock 
which performs, due to fluctuations, a biased random 
walk with velocity Vs (Fig. ^). If the incoming current 
j_ exceeds the outgoing current j+ , the shock velocity is 
negative. This is analogous to the back-moving shock of 
a traffic jam for sufficiently high incoming traffic fiow. 
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ti 



(a) 



XO + Vc{t2 - ^l) 
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P+ 



P-- 



(b) 



FIG. 3. (a) Diffusive spreading of a density perturbation 
in the stationary state at two times t2 > t\. The collective 
velocity describes the motion of the center of mass of the 
perturbation, (b) Motion of a shock. To the left of the domain 
wall particles are distributed homogeneously with an average 
density p_. To the right of the domain wall the background 
density is p+ > p^ 

To get an intuitive understanding how these veloci- 
ties determine the stationary phase diagram of a driven 
system coupled to boundary reservoirs let us assume 
P± = Pr,l- We consider first p+ > p* where p* is the 
density where the current takes its maximal value j* . 
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To make the argument more transparent we also assume 
that particles hop only to the right and that initially the 
lattice is empty. Because of ergodicity the stationary 
distribution is independent of the initial state and hence 
this assumption involves no loss of generality. Consider 
now the time evolution of the averaged density profile in 
a large system, starting from the empty lattice. We start 
the discussion by assuming the left boundary density to 
be very low. 

(i) As time proceeds, particles from the left reservoir 
will enter the system and (possibly after some distance 
describing the non-universal boundary layer) create a re- 
gion of constant density p_. This region decays to the 
right to zero, because after a finite time the rightmost 
particles will have traveled only a finite distance. Eventu- 
ally however, after a time which is of the order of system 
size, the rightmost particles will hit the right boundary 
with the reservoir of density p+. This reservoirs makes 
it more difficult for particles to travel further and hence 
creates a little traflSc jam. The result is a shock pro- 
file, with shock densities p-_ on the left and p+ to the 
right resp. like in Fig. ||b. [For the sake of the argu- 
ment one could have chosen such an initial state. The 
reason for choosing an empty initial state becomes clear 
below.] The decisive question is now how this shock pro- 
file evolves in time. According to the shock velocity 
under the circumstances described here is positive, sim- 
ply because the incoming current of particles j_ is less 
than the outgoing current for sufficiently small left 
boundary density. Hence, even though a shock forms 
by fluctuations, it has an average drift towards the right 
boundary. Hence the system remains in the low density 
(LD) regime with bulk density p = p- < p* ■ 

(ii) The situation changes when the left boundary den- 
sity takes a value such that the incoming current j_ 
equals the outgoing current In this case the shock 
velocity vanishes and the shock performs an unbiased 
random walk over the lattice. Hence the density profile 
may be regarded as being composed of two stationary do- 
mains with densities p_ and separated by a "domain 
wall" which is the sharp transitional region of the shock. 
Since the shock motion is unbiased, the stationary proba- 
bility of finding the shock is constant in space. This leads 
to an equal superposition of shock profiles and hence to 
a linearly increasing stationary density profile. In anal- 
ogy to first order equilibrium phase transitions where a 
domain wall separates regions of coexisting equilibrium 
regimes, we call the line defined by j+ = j_ and p_ < p*, 
/9+ > p*, a first order phase transition line. 

(iii) This line marks the transition to a high-density 
phase (HD) with bulk density p = p^ > p* since for even 
higher left boundary density the incoming current into 
the shock exceeds the outgoing current, and according 
to (H) the shock moves to the left. This leaves the bulk 
in the high-density regime determined by the coupling to 
the right boundary reservoir. At the first-order transition 
line the stationary bulk density is discontinuous, it jumps 
from p_ to p+. 



Next we consider the case of low right boundary den- 
sity p+ < p* ■ For definiteness we choose p+ = 0. The 
essential part of the following discussion, the explanation 
of the occurrence of a continuous phase transition, is un- 
affected by this choice. Again we start the discussion 
with the empty lattice. 

(i) As argued above, after some initial time the system 
will have filled up to a bulk density p = p-. Since p+ = 0, 
particles hitting the right boundary can leave the system 
without creating a shock. As a result, the system is in 
the low-density phase. [For p+ > a shock could form, 
but would have positive velocity under the circumstances 
considered here. Hence also in this case the system is in 
the low-density phase.] Now we examine the seemingly 
trivial reason why an increase in the left boundary density 
leads to an increase in the bulk density. Above we simply 
claimed this to be true on the basis of plausibility. Here 
we support this claim with an argument that becomes 
important below. Suppose we create a little perturba- 
tion at the left boundary by injecting an extra particle 
(on top of those particles that are injected anyway from 
the reservoir). This creates a perturbation which, by 
definition of the collective velocity, travels with Vc > 
into the bulk and leads to a local increase of the density, 
moving away from the boundary and spreading out as 
time goes on. The collective velocity is positive, since by 
assumption we have a background density p- < p* and 
hence the current as a function of the density has positive 
slope. The point is that maintaining such a perturbation 
(which corresponds to increasing the left reservoir den- 
sity permanently) leads to a permanent additional flow 
of particles into the bulk and hence to the anticipated 
increase of the stationary bulk density. 

(ii) It is now clear that this argument holds only as long 
as p- < p* ■ Assume now p- = p* ■ Following the reason- 
ing above this results in a maximal-current bulk density 
p* . However, if the left reservoir density increases beyond 
the maximal current density, the collective velocity be- 
comes negative. No extra particles flow into the system, 
which therefore remains in its bulk at the maximal cur- 
rent density p*. The system is now in the maximal cur- 
rent phase for all p- > p* and p+ < p*. This transition 
is continuous, the bulk density approaches p* smoothly 
from below. Intuitively this phenomenon may be un- 
derstood as "overfeeding" |^]. The injected particles act 
as blockages for further incoming particles, leading to a 
back-moving traffic jam at the origin. This increased den- 
sity at the origin blocks further injection attempts and 
prevents the actual increase of the current. 

Finally, using similar arguments, one can show that the 
transition from the high-density phase to the maximal 
current phase is also continuous. 

To summarize, the theory predicts a first-order tran- 
sition along the line defined by j+ = j_ and p_ < p*, 
> p* where the stationary bulk density jumps from p- 
(low density phase LD) to p+ (high density phase HD). 
On the phase transition line the stationary density is lin- 
early increasing from yO_ to p+. From both phases there 
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is a continuous phase transition to the maximal current 
phase MC defined by p_ > p* and p+ < p* ■ In this phase 
the bulk density takes the maximal current value p = p* . 
These rules are encoded in an extremal principle for the 
current 



max j{p) 

- _ / pe[pn,Pi] 

' min j{p) 



for PL > PR 
for PL < PR. 



(10) 



derived in |L1| . It is worthwhile pointing out that this 
dynamical theory explains in mesoscopic terms the pre- 
dictions one would obtain by viewing the system from 
a coarse-grained hydrodynamic view-point 1^]. At the 
same time, verification of this scenario suggests the va- 
lidity of the hydrodynamic approach to the description 
of the large-scale dynamics of the particle system [Q by 
using the exact current-density relation. 



B. Coupling to boundary reservoirs 

To verify this scenario which correctly describes the 
exactly solvable ASEP with open boundaries one has to 
investigate our model in terms of the effective boundary 
densities p_ , p+ . Given two reservoir densities there is no 
general recipe how to eliminate the non-universal bound- 
ary effects that result in the effective boundary densi- 
ties /5_ , p+ . Hence these quantities are not easy to con- 
trol. Ideally, for purposes of theoretical investigation, 
one would like to construct an injection and absorption 
mechanism which leads to a constant density profile for 
a semi- infinite system so that pL = p- and pR ^ p+, 
i.e. the effective boundary densities are identical to the 
actual control parameters of the model. 

Here we choose an injection mechanism where the par- 
ticles on the lattice interact with the reservoir particles 
in the same way as among each other. We define two 
injection rates from the reservoir at the left boundary: 

• injection at site 1 if site 2 is occupied 
|0A \AA with rate ai 

• injection at site 1 if site 2 is empty 
1 00 |A0 with rate 02 

and two new hopping rates at the right boundary: 

• hopping from site — 1 to site A^ 
yl0| 0A| with rate /3i 

• hopping out from site A^ (absorption) 
A] ^ 0| with rate /32. 

These four hopping rates are those that would be af- 
fected by the interaction with particles in the reservoirs. 
The rates ai {j3i resp.) have now to be determined as 
functions of the left (right) reservoir density. This can 
be illustrated e.g. for the injection process with rate 
ai. We imagine the reservoir to include a site of the 



chain. The injection rate into the first site is defined by 
the (stationary) average occupation of the imaginary 
site 0, but with the condition that the first site is empty 
and the second site is occupied. Considering the zeroth, 
the first and the second site as three neighboring sites of 
an infinite chain, this conditional probability can be ex- 
pressed readily as correlations in the stationary state of 
an infinite chain. Thus we find ai = q{ 101 )/( 01 ) where 
expectation values like (101) = ( 71^(1 — ni+i)nj+2 ) are 
calculated in the thermodynamic limit of the distribution 
(^ with density p = pL- The case of a2 is entirely similar 
and one finds a2 = r(100)/(00). 

For the calculation of the right boundary rates we note 
that the probability of jumping from the site A — 1 to 
the site A^ is affected by the average occupation of the 
imaginary reservoir site A" -I- 1. With the jump condition 
that site A^ — 1 is occupied and site A^ is empty one finds 
(3i = (r( 100) + q{ 101))/( 10). The case of (32 is similar 
but one has to take into account the conditional proba- 
bility of the occupation of two imaginary reservoir sites 
A + 1 anAN+2. This gives /32 = (r( 100 ) -hg( 101 ))/( 1 ). 
One has to determine these correlations in an infinite 
chain with density pR. These rates can be expressed as 
a function of the density through the form of the current 
and we find: 
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For PL = 1 we use the limiting values ai = q, a2 = r, 
and for pR = we use fii = (^2 = r. Together with the 
bulk hopping rates r, q the dynamics of the model is now 
completely defined. 

Before discussing the full phase diagram we consider 
the case of equal boundary densities p — Pr — Pl- Some- 
what surprisingly, it turns out that we can actually ob- 
tain the full stationary distribution of the process: 



p*(M) = 



i-P 



EN-1 

= 2 "'(Ai-l)"i+"" 
(16) 



with the eigenvalue Ai of the transfer matrix of the one- 
dimensional Ising model and the "fugacity" z = e~^^ 
(Appendix A). Stationarity of this distribution can be 
proved by writing the time evolution operator of the 
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process in a quantum Hamiltonian formalism (see ap- 
pendix B). Like in the periodic case, the exact station- 
ary non-equilibrium distribution of the open system with 
equal boundary densities is the equilibrium distribution 
of an Ising chain of length N, but with boundary fields 
g{ni + njv) rather than with the Ising coupling jintii. 

Moreover, with the transfer matrix formulation of the 
Ising model it is straightforward to show that the den- 
sity profile is constant. The theoretical scenario de- 
scribed above then suggests the identification p_ = 
for PL < p* and p+ — pn for p^ > p*. For p^ > p* or 
Pfi < p* respectively constant profiles are not stable with 
respect to fluctuations and the relationship between the 
reservoir densities and the effective boundary densities 
may be more subtle (see below). 



IV. PHASE DIAGRAM FOR REPULSIVE 
INTERACTION 

Having defined the model and given the theoretical 
background we are now set to investigate the phase di- 
agram of the system. We consider first the traffic flow 
scenario (repulsive interaction). 



A. Mean field 

Unfortunately the exact stationary distribution of the 
open system for different boundary densities does not 
have a simple form. However, also the density in the 
open system satisfles a continuity equation of the form 
(H), with the boundary currents 

jo = ai( (1 - ni)n2 ) + (1 - "•i)(l - "-2) ) (17) 
jAT-i = /3i(riAr_i(l - n^v)) (18) 
jN=02{nN) (19) 

Using also the expression (^) for the bulk current one 
can obtain the stationary density profile in a mean field 
approximation by neglecting the correlations in the ex- 
pectation values. Stationarity implies equal current ev- 
erywhere in the lattice, (ji) = j. This gives rise to the 
bulk recursion relation 



Pi = 



(1 - Pi+i) {qPi+2 + r{l - p.1+2)) 



(20) 



for the density profile. The boundary values of the re- 
cursion are determined by the boundary currents (|l7| ) - 

Because of the interaction with the boundary there 
are two possibilities to perform a mean-field analysis. 
The most straightforward choice would be to choose the 
boundaries rates ( pi] ) - (14) and analyze the mean field 
phase diagram in terms of the reservoir densities pr^l- 
For a given choice of current one can draw a density pro- 
file and read off the various phases. However, it turns 



out that within such an approximation scheme a constant 
density profile cannot be achieved even if p_R = pl and 
that the boundary densities are not equal to the reser- 
voir densities. As a result, the mean field phase diagram 
obtained in this way differs considerably from the theo- 
retical expectation. Within this approximation the MC 
phase disappears for r / q smaller than « 0.6. 

From a theoretical point of view it is more natural to 
neglect the correlations even in the expressions of the 
boundary rates in terms of the Ising expectation val- 
ues. This gives rise to simple expressions of the boundary 
rates in terms of the reservoir densities. Using the bound- 
ary conditions po = pl and pn+i = Pn+2 = Pb. extends 
the validity of the recursion ( pO| ) to the range < i < N . 
In this way pN can be expressed as a function of pn and 
j and this was the reason of choosing the direction of the 
recursion from the right to the left. It is easy to see that 
for a semi-infinite system there are the constant solutions 
Pi = pi^ and Pi — pr resp. with the mean field current 



f '^{p) = p{l-p) {qp + r{l-p)) 



(21) 



Like in the exact solution this is also a solution of the 
finite system with equal boundary densities. Therefore 
we shall discuss this mean-field approach in some more 
detail. 

For general boundary densities we did not find a so- 
lution of the recursion (^o|) in closed form, but it is eas- 
ily solved numerically on a computer. It is sufficient to 
choose a density pji and a current j and apply the re- 
cursion relation for drawing the density profile. If any 
of Pi's is smaller than or greater than 1 then there is 
no solution corresponding to these values of pR and j at 
this system size since any physical solution has to satisfy 
< Pi < 1 for all 1 < z < A^. For fixed pr,pl and 
length N this requirement fixes the current j. In this 
way one can map out the whole phase diagram, finding 
the corresponding current for all values of pL and pR. 
For sufficiently large N the size dependence is negligibly 
small. 

We have chosen r = 1.0 and q = 0.1. The low- (LD) 
and the high-density (HD) phases resp. are easy to dis- 
tinguish as the bulk density is equal to one of the end den- 
sities. This allows us to identify p_ = p^ and p-^- — pR. 
At the other end oscillatory behavior is observed. This 
does not happen in the usual ASEP with open boundaries 
but has been observed in an ASEP with particles 
covering more than one lattice site The maximal 
mean field current accessible at given parameters r and 
q and the corresponding bulk density define the maximal 
current phase MC. The location of these phases in the 
parameter space differs from those obtained from theo- 
retical scenario reviewed above (Fig. ^ by using the exact 
value (|^) of the current. Further analysis shows that the 
discrepancy increases with increasing repulsion. 
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Repulsive case 
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FIG. 4. Mean field phase diagram (broken curves) and the- 
oretical phase transition lines (full curves) for repulsive inter- 
action. 



B. Monte Carlo results 

In order to check the phase diagram based on the 
theory of boundary induced phase transition described 
above we have performed Monte Carlo simulations of the 
model. 

The boundary rates are determined as explained above 
to simulate given left and right densities. During one MC 
step one of TV -I- 1 sites is chosen randomly {N is the size 
of the system and the +1 is for the jumping into the 
system) and if there is a particle then it can hop with 
a probability given by the hopping rates. The initial 
configuration was an empty lattice. The required time 
to reach the stationary state for given system size can be 
determined by investigating the time dependence of the 
current and the bulk density (the density in the middle of 
the system) . The order parameter which is the stationary 
value of the bulk density obtained as an average over 10^ 
MC steps is the best indicator of the phase transition 
lines. At the first order phase transition between the 
high and the low density phases it has a pronounced jump 
from the left to the right density. A system size of 1000 is 
sufficiently large to localize the phase transition line. The 
transition into the maximal current phase is continuous, 
only the derivative of the order parameter has a (jump) 
discontinuity. In order to localize this phase transition 
line one needs larger systems of size 5000. 

In a given phase the simulated bulk density depends on 
the boundary densities the same way as it is predicted by 
the theory, namely it equals to the left (right) density in 
the LD (HD) phase, and it is a well defined constant value 
in the MC phase. At the phase boundary the crossover 
between two kinds of behavior of the bulk density can 
be localized within an error shown the figures (). The 
Monte Carlo results are seen to be in agreement with 



the phase diagram derived from the theory of boundary- 
induced phase transitions (Fig. The mean-field phase 
transition lines are significantly outside the numerical er- 
ror bars except on part of the transition line between LD 
and MC phase. 
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FIG. 5. Monte Carlo data and theoretical phase transition 
lines (full curves) for repulsive interaction. The numerical 
error bars are indicated by full circles. 



V. ATTRACTIVE INTERACTION 

Attractive (ferromagnetic) interaction leads to bulk 
particle domains as in the one-dimensional Ising model. 
Dynamically this is qualitatively understandable since 
because of g > r particles tend to form clusters. One ex- 
pects a shift of the maximal-current density p* to higher 
densities. This can indeed be shown by calculating the 
derivative of the exact current (^. The full current- 
density relation for q = 1, r = 0.1 is shown in Fig. ^ 



Attractive case 
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FIG. 6. Exact stationary current j as a function of the 
density p for g = 1, r = 0.1 (repulsive interaction). 

Analysis of the model with open boundaries proceeds 
in the same way as in the repulsive case. With r = 0.1 
and q = l.O the mean field approximation yields the three 
phases discussed above, with the phase transition lines 
differing substantially from those expected theoretically 
if one identifies p_ — and pj^ = pj^ as suggested by 
the density profiles in the low and high density phases 
respectively. 

As additional feature the mean field theory predicts 
a further phase transition located at the large pi^ and 
small PR corner of the phase diagram (Fig. |7|). In this 
"fourth phase" the bulk density is larger than the density 
corresponding to the maximal current. The area of this 
phase increases with increasing attraction, and the area 
of MC phase decreases but does not disappear. 
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FIG. 7. Mean field phase diagram (broken curves) and the- 
oretical phase transition lines (full curves) for repulsive inter- 
action. The narrow area above the HD phase is the mean field 
MC phase. The fourth phase obtained from mean field covers 
the large area in the upper right part of the phase diagram. 



terms of which the theory is formulated). While for re- 
pulsive interaction and the choice of coupling mechanism 
used here these two quantities may be identified, the con- 
nection is apparently more complicated for attractive in- 
teraction if the (real) boundary densities become suffi- 
ciently high (left boundary) and low (right boundary) 
respectively. An adequate explanation of this boundary 
phenomenon which leads to a reentrance transition into 
the high-density phase is not available. In order to rule 
out the possibility of a finite-size effect we performed sim- 
ulations with different system sizes. Using N=500, 1000 
and 5000 resp. we found no indication that the reen- 
trance transition would disappear for large enough sys- 
tems. 

Attractive case 



a: 
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FIG. 8. Monte Carlo data and theoretical phase transition 
lines (full curves) for repulsive interaction. The numerical 
error bars are indicated by full circles. The reentrance phase 
transition to the fourth (HD) phase (upper right corner of the 
phase diagram) is marked by the broken curve interpolating 
between the Monte Carlo points. 

VI. CONCLUSIONS 




The phase diagram obtained from Monte-Carlo simu- 
lation has the phase transitions expected from the theo- 
retical prediction, but also shows that the fourth phase 
exists in a small neighborhood of the {p^ — 0, p^ > p*) 
line and hence is not an artefact of the mean-field ap- 
proximation. The location of the phase transition lines 
is rather different compared to mean field (Fig. ||). The 
transition to this phase from the LD phase is first order 
and from the maximal current phase it is second order. In 
this sense the fourth phase is like the high-density phase, 
however, the bulk density differs from both the boundary 
densities pr^l- Within the theory of boundary-induced 
phase transitions this high density phase can be under- 
stood by recalling the non-universal relationship between 
real boundary density and effective boundary density (in 



Our analysis of the totally asymmetric exclusion model 
with next-nearest-neighbor interaction consists of two 
parts. First we considered the periodic system in or- 
der study the bulk properties of the steady state. With 
a view on traffic flow modelling we remark that the ob- 
servations (i) - (iv) listed above are consistent with the 
results obtained here. The current-density relation be- 
comes asymmetric (Fig. |^) in a way which is closer to 
real traffic data as the symmetric relation j — p{l ~ p) 
for the ASEP with r = q = 1. There is no discontinu- 
ity in the derivative at the maximal-current density p*, 
in agreement with the arguments given above. As more 
far-reaching conclusion we note that the exact result (Q) 
sheds light on the so far somewhat unclear relationship 
between the updating mechanism and the occurrence of 



9 



correlations. In order to obtain correlations in traffic flow 
models it is not essential to use a discrete-time updating 
mechanisms. Since other discrete-time models have un- 
correlated stationary states it appears that correla- 
tions have their physical origin in speed-reduction rather 
than in the nature of the updating scheme. 

In the second part we investigated the steady-state se- 
lection in the open system. It turns out that the theoret- 
ical scenario based on the interplay of shocks and over- 
feeding correctly describes the phase diagram in terms of 
effective boundary densities. For repulsive interactions 
this strengthens the case for a maximal current phase 
in real traffic. For attractive interaction there is, how- 
ever, a surprising reentrance phase transition to the high 
density phase which originates in the so far poorly under- 
stood relationship between actual and effective boundary 
densities. Apparently these two are not always monotic 
functions of each other as one would naively expect. For 
a deeper understanding the next step to be done is the 
analysis of the density profiles close the reentrance phase 
transition lines in order analyze whether universal prop- 
erties of the usual transition lines ||,^,|l^ can be ob- 
served or not. 

T.A. thanks for partial support by the Hungarian 
Academy of Sciences (Grant OTKA T 029792). G.M.S. 
thanks K. Klauck, V. Popkov, L. Santen and A. Schad- 
schneider for illuminating discussions on traffic flow. 



APPENDIX A: TRANSFER MATRIX FOR THE 
ONE-DIMENSIONAL ISING MODEL 

The energy of the one-dimensional Ising model with 
the classical spin variables Si = ±1 can be written 
in term of the variables = (1 — Si)/2 in the form 
E = JJ2i nitii^i. In this interpretation the Ising model 
is a classical lattice gas of M = hard-core particles 

with nearest neighbor interaction. Using standard tech- 
niques p3] one can express the grand canonical partition 



function Z = 



con fig 



pE inverse temperature 
(3—1/ kT in terms of the eigenvalues of the transfer ma- 
trix 



1 ^/I 



z 



(A.l) 



where, for our model, e~^'^ — q/r and we define z = 
e~^^. In the spin interpretation of the model h plays the 
role of a magnetic field. 

The eigenvalues of the transfer matrix are 



X,^, = \{l + lz)±^\{l + lzY + z{l-l) (A.2) 

where we choose Ai to the positive sign. For a sys- 
tem with N sites and periodic boundaries one has Z — 
TiT^ = + X2 ■ The equilibrium distribution of the 



Ising model is the stationary distribution of the parti- 
cle hopping model. This does not mean that the parti- 
cle hopping model reaches thermal equilibrium at long 
times, since the stationary distribution does not satisfy 
detailed balance with respect to the dynamics of the 
model. The non-equilibrium nature of the steady state 
results in a non- vanishing stationary particle current. We 
stress that completely different dynamical models may 
have the same stationary distribution, see e.g. the list 
of models given in ||2l[| . Another non-equilibrium ex- 
ample with the same Ising distribution is the discrete- 
time ASEP with parallel update Even though the 
distribution is the same and hence all particle correla- 
tions are the same, the stationary current of this model 
IIitI , ^ (which is defined via the dynamics by the conti- 
nuity equation) is different from the current (^) in our 
model. 

In order to calculate expectation values in the transfer 
matrix formalism we define the diagonal matrix 




1 



(A.3) 



The density is then given hy p — Tr{nT )/Zn and 
a two-point correlation function is given by ( Ufen; ) = 
Tr(nT''~'nT^~''+'). Higher order correlators are calcu- 
lated analogously. By diagonalizing T one obtains the 
expression 



1-Ai 
A2 — Ai 



(A.4) 



for the particle density in terms of the eigenvalues of T. 
This relation may be used to express the current (||) as 
a function of the particle density. In the thermodynamic 
limit N —I- 00 one obtains (^. 

We note that expectation values for the distribution 
( p^ ) are calculated with the same transfer matrix as in 
the periodic case. However, instead of taking a trace, one 
calculates a scalar product with suitably chosen vectors 
which are determined by the boundary fields. 



APPENDIX B: PROOF OF STATIONARITY FOR 
EQUAL BOUNDARY DENSITIES 

The proof of stationarity of the distribution (|^) for 
the periodic system is given in Ref. An alternative, 
constructive proof can be obtained by using translational 
invariance and taking the distance between neighboring 
particles as stochastic variables. In this way the particle 
hopping process turns into a zero-range process for 
which the stationary distribution is known. Reexpressing 
the stationary zero-range distribution in terms of particle 
occupation numbers yields (^. 

Here we prove stationarity of the Ising distribution ( p^ ) 
for the open system coupled to reservoirs of equal density. 
We use a convenient standard approach for the stochas- 
tic description of classical interacting particle systems. 
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known as "quantum Hamiltonian formalism" Q. The 
basic idea is to formulate the generator of the Markov 
process in terms of a many-body quantum operator. For 
processes without exclusion one obtains in this way a 
Fock space representation of the generator in terms of 
bosonic creation and annilation operators |37-39|. For 
exclusion processes with at most one particle per site the 
same strategy yields an operator expressed in terms of 
Pauli-spin matrices 40 - 4^] . 

We define the exclusion process with state space X = 
{0, 1}^ and transition rates Wn^n' from state n to n' in 
terms of a master equation 

j^P{mt)^ [w„.^„P(ri';t)-w„^„,P(n;i)] (B.l) 



for the probability P(n; t) of finding, at time t, a con- 
figuration n of particles on a lattice of sites. Here 
n = {ni, 712, . . . , nAf} where Ui =0,1 are the integer- 
valued particle occupation numbers at site i. The master 
equation is a linear first order DGL in the time-variable 
and therefore it is natural to write it in a vector notation 
with the probabilities P{n] t) as vector components. We 
represent each of the possible particle configurations n by 
a column vector | n ) which form a basis of a vector space 
X — (C^)^^. The transposed vectors {n \ form a basis 
of the dual space and we define the usual scalar prod- 
uct {n\ri!) = 5n,n'- The probability distribution is now 
represented by a state vector | P{t) ) = X^nGX ^(b; t)\ll) 
and one can write the master equation in the form 



dt 



Pin;t) 



[n\H\P{t)) 



where the off-diagonal matrix elements of H are the (neg- 
ative) transition rates between states and the diagonal 
entries are the inverse of the exponentially distributed 
life times of the states. In formal analogy to the quan- 
tum mechanical Schrodinger equation we shall refer to H 
as quantum Hamiltonian. A state at time t' = + t is 
given in terms of an initial state at time to by 



\pito+t): 



-Ht\ 



(B.3) 



We stress that the physicists notion "quantum Hamil- 
tonian" for the matrix H is somewhat misleading in so far 
as H is, in fact, the generator of the Markov semigroup 
of the process, rather than the Hamiltonian of an actual 
quantum system. This by now well-established notion 
has its origin in the fact that for various stochastic pro- 
cesses the generator H is identical to the quantum Hamil- 
tonian of some well-known spin system. In this context 
we would also like to point out that quantum mechanical 
expectation values {A) = ( ^I* 4* ) for an observable A 
are calculated in a different way than probabilistic ex- 
pectation values for a function F{n) of the stochastic 
variables 77. In the quantum Hamiltonian formalism one 
writes {F) = J2nex t) ^ {s\F\ P{t) ) with the 



matrix F = X^nex Pilii)\ll){ll \ t;he summation vec- 
tor ( s I = X^nex I which performs the average over all 
possible final states of the stochastic time evolution. 

For our considerations the expectation value pk{t) = 
{s\nk \ P{t) ) for the density at site k is of special inter- 
est. It is given by the projection operator nk which has 
value 1 if there is a particle at site k and otherwise. 
TO-point density correlations are then given by the ex- 
pression ( s |nfcj . . . rik,^ I P(t) ). In this paper we study 
only stationary expectation values. For the formal de- 
scription of a stationary probability distribution we use 



the transposed summation vector | s) = XnGX I—)- ^ 
general stationary measure P*{n) may then be written 
in vector notation in the form | P* ) = e~^^^-^ s)/Zi\[ 
with the configuration-dependent "energy" matrix E{n) 
and the "partition function" = (s \e~^^^-^ \ s) which 
normalizes the measure | P* ) = e^^^^^-^ \s). Notice that 
in vector notation the expression E{n) represents a di- 
agonal matrix with the energies as diagonal elements. 
Below we shall also use the invertible diagonal matrix 
P* — e~^^'-^ which has the (unnormalized) stationary 
probabilities P*{n) as diagonal elements. 

To obtain the quantum Hamiltonian for the time evo- 
lution of the interacting ASEP with open boundaries we 
note that one can represent any two-state particle system 
as a spin system by identifying a particle (vacancy) on 
site k with a spin-down (up) state on this site. This al- 
lows for a representation of H in terms of Pauli matrices 
where nk = (1 — crfe)/2 projects on states with a particle 
on site k and Vk — 1 — nfc is the projector on vacan- 
cies. The off-diagonal matrices = (tr^ ± ia^)/2 create 
(B.2) (s^) and annihilate (s^) particles. We stress that in the 
present context the "spins" are just convenient labels for 
particle occupancies. Using this pseudospin formalism 
one finds 



N-2 



^ = X! i^kVk+1 - s^Sj.^ J(rnfe+2 + qVk+2) 

k=l 

+ {1 - ni - s^){ain2 + a2V2) (B.4) 
+Pi{nN-iVN - sjj^is]^) + P2inN - s%). 

Within this formalism proof of stationarity is now a 
straightforward calculation, using simple expressions for 
the jumping rates at the ends: 



ai = 

a2 = 

/32 = 



qz 



Ai(Ai-l) 

rz 

Ai(Ai-l) 

Ai-l 

Ai' 



(B.5) 



and assuming that the bo unda ry densities are equal at 
the two end. According to ( B.3 ) stationarity is equivalent 
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to the eigenvector relation H\P*) = which in turn is 
equivalent to 

{P*)-^H\P*) = {P*)-^HP*\s) ^ 0. (B.6) 

The diagonal similarity transformation of H with P* 
leads to a sum of transition matrices which act non- 
trivially on at most four neighboring sites. To calcu- 
late the action of the non-diagonal parts on | s ) we use 
s'^l s) = Vk \ s) and s) — nk \ s). This leaves only di- 
agonal terms the sum of which vanishes identically. This 
proves stationarity of the measure (|l6|). 
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